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1. Introduction 


At present in computational gas dynamics [1] there are well- /3 * 
developed methods based on isolation of discontinuities, making 
use of calculation grids adapted to the boundary. Such approach 
[2-3] affords high accuracy, although the computerized construc- 
tion of a difference grid involves major software difficulties, 
especially in the case of a moving boundary and other nonsta- 
tionary processes. 

The tasks of contemporary engineering gas dynamics are dis- 
tinguished by great diversity of geometries and boundary condi- 
tions, and even an approximate solution often enables an assess- 
ment as to the prospects of a gas dynamics apparatus under 
development. There also exist a multiplicity of problems in 
which the gas flow itself plays a subordinate role in relation 
to the other dominant physical processes for which only extremely 
crude models are available (physical gas dynamics). In these 
applications, it seems advisable to expect the following demands 
of the computer algorithm: 

1. A unified approach to the solving of different types of 
problem, with easy enlargement of the physical content. 


* 
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2. Easy operation of the program, not requiring a highly 
qualified computer mathematics user. 

3. Satisfactory practical precision at relatively short 
computer run time. 

It is difficult to satisfy requirements 1-3 with traditional 
numerical methods. The methods best suited to this are a start-to- 
finish calculation [1] , which has been a stimulus to the develop- 
ment of the latter. The use of boundary-nonadapted regular grids 
and explicit schemes enables a simplification of the algorithm, 
an automated running of the program, and consequently an appre- 
ciable enlargement of the range of users. 

Various authors have elaborated this method for the two- 
dimensional case [4-9], The chief shortcomings of the methods 
described in [4-6] are the nonconservativeness , the nonhomo- 
geneity and the cumbersome formulas in the boundary compartments 
of the grid, which do not allow a generalization of these methods 
to the case of an arbitrarily moving, curvilinear boundary in a 
complex three-dimensional configuration. This complicated algo- 
rithmic procedure has another shortcoming — the boundary condi- 
tions are difficult to realize on nonuniform rectangular grids. 

The methods advocated in [5,6] are based on direct approximation 
of the derivatives at the node next to the boundary, which 
results in a large error at nodes with a heavily nonuniform 
pattern, complicates the software realization, and inadmissibly 
limits the time interval of explicit layouts. These shortcomings /4 
led the author [5] to the conclusion that nonadapted grids have 
few prospective applications. 

However, the above does not apply to the scheme [7-9] derived 
on nonadapted grids by the integro-interpolation method (IIM) 

[10] . Use of the IIM enabled simple difference formulas for the 
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two-dimensional case, assuring conservativeness of the start-to- 
finish calculation of the flow parameters in the interior and 
boundary compartments even for moving boundaries of complicated 
shape. Joining small compartments to the neighboring ones pro- 
vides adequate stability of the difference scheme. 

The present work generalizes the approach developed in 
[7-9] to the case of three-dimensional solids of complicated 
shape. 

2. The Numerical Method 


Let us consider the equations describing multidimensional 
flow of a nonviscous gas in Eulerian variables: 


+ o/< v <§ WJ .-0 
+ 1) tv^W.W/* jiad P - Q 



+ c/iv- WJ + oliv (P-ty) -zQ 

P= P(s,£) 


( 1 ) 


or the equation for the total energy E can be replaced by the 
equation for the internal energy e: 


i — . 

4 - c1i*{g£W) + P. w=o 

where W = (W ,W ,W ) is the velocity, P the pressure, p the 
density. System (1) is solved in the region G with boundary 

0 I? 

r = r' 1 ' (Fig. 1), where T' 1 ' can be of three kinds: 1) the entry 

£ 

boundary , where conditions of the first kind are usually 
set; 2) the exit boundary , where we shall set the condition 
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of vanishing of the normal derivatives of the gas dynamic quan- /5 

Z Z — 

tities; 3) the solid wall or planes of symmetry r , where we 

have the condition of no flow-through: 

I 

- f>) - o j w 

where n is the normal to the wall. 

The surface of the solid is given by the function F: 


F(x,y,z) > 0, point (x,y,z) lies outside the solid 
F(x,y,z) = 0, point (x,y,z) is on the solid (3) 

F(x,y,z) < 0, point (x,y,z) lies inside the solid. 


z z z 

The boundaries T, , T, „ and r are situated in the coordi- 

DX O^IX O 

nate planes. We shall complement the region G to produce a 
coordinate rectangle by introducing the region G T , corresponding 


to the solid. 

rectangular grid, producing three types of compartment 


We cover the rectangle G + G^ with a nonuniform 


a) a full compartment, entirely situated in the gas (inside 
region G) 


b) a partial compartment, intersected by the boundary of 
the solid (in the calculations of the field of flow only that 
portion occupied by gas is involved) 

c) an empty compartment, entirely situated in the solid 
(inside region G T - hereafter omitted from the calculations, and 
such compartments are totally absent from the specially organized 
computer memory files. 



Furthermore , 
along r£„ x , 


we introduce a layer of imaginary compartments 
After this, the grid is complemented to a 
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rectangle by the region G^ 0/7 , the compartments of which are 
equivalent to empty ones. The center of each compartment corre- 
sponds to the index i = (i^,^,^). Hereafter, we shall omit 
from the variables the indexes of directions on which no differ- 
ence operations are performed. For the compartment of type a,b 
(Fig. 2), we introduce the following geometrical parameters: 

S+ = - surfaces of the side faces of the compartment i 

open to passage of gas in the direction i = 1,2,3; - the 

volume of the compartment occupied by gas. The proposed tech- /6 

nique can be developed on the basis of any given explicit diver- 
gent (flow) method of start-to-f inish calculation. The results 
of the two-dimensional calculations [7-9] revealed that the FLIC 
method [11] is a convenient reference point, in which the values 
of all the gas dynamic variables <p = [p,pW,pE] are adjusted to 
the centers of the compartments. 

This method is based on consecutive calculation of the 
physical processes - transition from the n-th to the (n+l)-th 
time layer is done in a two-stage (in the case of moving bound- 
aries, three-stage) scheme with a summary approximation [10] . 

In the first stage, only the action of pressure forces 
with transport processes held in abeyance is considered. To 
smooth out pulsations behind the shock wave front, it is possible 
to introduce an artificial viscosity. We shall assume such 
viscosity by analogy with the second physical [viscosity?] [9] : 

in contrast with [4,5, 7,8] this is added to the gas pressure [3]: 

p'= P + ? 

which simplifies the algorithm and improves the equalization of 
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pulsations in the multidimensional case (hereafter we shall 
write the combined pressure without the prime) . Let us apply 
the procedure described in [7] to a compartment or test volume. 
We integrate system (1) throughout the volume of the compart- 
ment i and replace the time derivative with its difference 
analog : 


(<^ W ) - (<S WJ + ^ g r acl P cJV =0 


v W JV*0 


(4) 


The index 0 refers to the value of the variables in the n-th 
layer, 1 after the first stage, 2 after the second. Using the 
familiar formulas of vector analysis, we change the volume inte- 
grals in (4) into surface integrals: 

4 (§v7) = ’(swi +-r ^P-ct£ 

\ S £) = # (§f) -tr-Pj (5) 

o 


On the side faces of the compartment open to the passage of gas /7 

all the quantities shall be approximated by the half-sum of the 
values in the corresponding compartments. It follows from (2) 
that (W’dJa) = 0 on the solid wall: therefore, the wall in (5) 
provides a zero contribution to the energy integral. Provided 
that the spacing of the grid is much less than the radius of 
curvature of the solid, we have from (2) : 

(6) 

which lets us regard the pressure on the wall in (5) as equal to 
the pressure at the center of the compartment. 
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As a result, we obtain approximation formulas for the first 
stage : 



(‘w“- -W + 

(*£ - ’£) §V = - p.t . ± l s ,m *W.i - 
- S'., • w. t ) 

*W= 0.S'(°v7 + *W) 
v/.l »o.r- <>«•%■ W“) 


i 

I 


In the second stage, convective transport is factored in. The 
procedure of obtaining the difference formulas is similar to the 
first stage procedure. In this case we use a "scheme with wind- 
ward differences". This scheme conveys rather precisely the 
characteristic physical features of the flow, i.e., it possesses 
the attribute of transportiveness. For the second stage, we get 
the following difference formulas: 


a- 't = os s’" a' w.j 

* (*wn 

* 2 


For the first-order explicit method described, the Courant con- /8 

dition is a criterion of stability (for more details cf . [12]). 

To avoid excessively severe limitations on the time interval, 
each small partial compartment (with volume less than 20-30% of 
a full compartment) is joined to its neighbor. 
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3. The Software Realization 


The program AEOL-3 is written in the language FORTRAN, 
making use of the conventions of the OLYMPUS system [13,14] and 
being an expanded version of the AEOL-2 application program 
package for the three-dimensional case. The structure of the 
program and the interplay between its modules conform to the 
general layout of the AEOL package. The programming does not 
employ any specific features of a particular FORTRAN version, 
and therefore can be used in virtually any computer of sufficient 
capacity. The package was developed in the BESM-6, but is run 
in the YeS-1045. The computer run time is relatively short: 
establishing a stationary flow past a complicated three-dimensional 
solid on a grid of 22,000 compartments (working storage around 
1 MGB) takes less than 3 hours of machine time of the YeS-1045. 

Only the construction of the modules CURVE and CUBE differ 
significantly from the two-dimensional version AEOL-2. These 
perform the computation of the geometrical parameters of the 
partial compartments. If the boundary is immobile, the sub- 
programs are called up only once at the start of the program. 

CURVE effects a scanning of the compartments of the grid, the 
geometrical situation in each partial compartment being solved 
independently, which also enables a uniformity of solving dif- 
ferent types of problem. The CURVE and CUBE modules employ a 
function for specification of the shape of the solid (3) , formu- 
lated as a subprogram - the function FCUR(x,y,z). The lengths 
of the edges of a partial compartment are determined by solving 
the equation FCUR(x,y,z) = 0 on the corresponding edge by the 
method of dividing the segment into halves. If the boundary of 
the solid passes exactly through one of the vertexes of the 
compartment, we shift the boundary a slight distance away from 
the vertex (usually less than 0.1%) for uniformity of the algo- 
rithm. 
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To simplify the finding of the area of a particular side 
face of a partial compartment, the line of intersection of the 
face with the boundary of the solid is approximated by a straight 
line segment. If (q^,q 2 ) an d are the pairs of lengths 

of opposite edges of the same face, it is easy to obtain a homo- 
geneous formula: 


S = 0.5*-(q t +qj • - v»i*(Q t Q t ) .*,<* (£ £] 

The subprogram CUBE calculates the volume of the portion of /9 
a compartment with a lesser number of vertexes, which is sub- 
tracted as needed from the volume of the corresponding full 
compartment. Through geometrical manipulations (rotation, 
mirror and central symmetry, parallel translation), the spatial 
configuration is reduced to one of the standards. Depending on 
the type of configuration, we calculate the volume by a particu- 
lar approximation formula, and if the situation is nonstandard, 
the program is automatically halted and an error notification 
is sent. (A nonstandard situation is easily corrected by modi- 
fying the grid, cf. Fig. 3f.) The version being used at a given 
time provides five standard configurations (Fig. 3a-e) , which 
enable analysis of practically all configurations. 

The described standardization of cases of intersection of 
a three-dimensional compartment by the solid substantially sim- 
plifies the program logic, reduces the run time, enables a uni- 
fication of the program complex and facilitates its mastery by 
the user having no special training. The user is only expected 
to specify the shape of the solid in the FCUR function program. 

On the foundation of a unified geometrical data base, two 
subsystems have been created in the AEOL-3 program complex: 
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a) calculation of the hydrodynamics, b) visualization 
on display screen of the shape of an intricate object by lighting 
up the lines of intersection of the surface described by the 
function FCUR(x,y,z) with the planes of the grid. This program 
appreciably simplifies the shape adjustment process. 

4. Test Computations 

To test the proposed procedure of start-to-f inish calcula- 
tion, we shall solve the problem of a supersonic gas flow past 
a particular standard solid on various grids and shall compare 
the results with standard tables. A sphere on a uniform Car- 
tesian rectangular grid is a solid of practically arbitrary 
shape (Fig. 3g) . Therefore, we examine the streaming of a gas 
flow with parameters = 2, y = 1.4, P oo = P ro = 1 around a unit 
sphere on grids 21 x 24 x 32 and 18 x 18 x 24 with constant 

spatial intervals h = h = 0.103 and 0.164, h = 0.0594 and 

x y z 

0.095, respectively. On the first grid, in the entire field of /IQ 
flow outside the zone of discontinuity the differences of the 
pressure and density from the tabulated data [15] did not exceed 
3% and 6%, while the drag was in agreement within 1.5%. 

Figure 4 shows graphs of the pressure distribution along 
radii making angles co = 0°, 60°, 90° with the direction of the 
flow; Fig. 5 shows the pressure distribution along the solid. 

The curves obtained from a calculation on a detailed grid are 
indicated by 1, those from a more coarse grid by 2. Figure 6 
and 7 show the corresponding density profiles. It is evident 
from the graphs that the agreement of the pressure and density 
fields on both grids is perfectly satisfactory. Comparison of 
the calculation results on the two grids demonstrates the con- 
vergence of the method as the grid becomes more fine. It should 
be pointed out that an accuracy better than 1% should not be 
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expected for realistic grids in our method, since the parameters 
within the grid compartment are assumed to be constant. In order 
to clarify the role of the mutual arrangement of the solid and 
the grid, the sphere flow problem was again solved on a grid of 
24 x 24 x 32, but with the center of the sphere shifted 0.5*h 

z 

along the flow. As compared to the former case, the difference 
in pressure field was less than 1%, but the difference in the 
drag force remained within 1.5%. This result is a dramatic illu- 
stration of the reality of the above-mentioned accuracy limit 
level of 1%. 


In problems of experimental design the foremost considera- 
tion is the distribution of pressure along an object and the 
aerodynamic coefficients, as well as the nature of their change 
as the shape of the object is varied. Analysis of the test 
calculations reveals that the developed technique, realized in 
the AEOL-3 program complex, is able to solve such problems on 
extremely coarse grids with tolerable precision. 

5. Calculation of the Flow Around a Three-Dimensional Object 
of Intricate Shape 


Let us consider the problem of a gas flow with parameters 
M oo = 3.5, p oo = p oo = l, y = 1.4 around an obtuse object. The 
general appearance of the object in Fig. 8a is shown by the 
lines of intersection of the object's surface with the coordi- 


n 3 +^0 pig. no c ( 
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■*- XM UJ. L/ UllNX >-< 
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the plane of symmetry and in the horizontal meridional section. 

About the periphery of the graph are also shown the grid plane 
markings, illustrating the distribution of the intervals. It 
should be noted that, despite the small dimensions of the recess 
in the meridional plane as compared to the grid interval (Fig. 8c), 
the method provides an excellent resolution of the secondary 
shock. The lift force at such object is negative, which was /II 
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found in the calculations and is entirely explained by the 
presence of the secondary shock, which increases the pressure 
in the upper half. 


The calculations were done on grids with different spatial 
intervals, the difference in the vertical and horizontal aero- 
dynamic coefficients not exceeding 2%. A standard analysis was 
done on a grid of 15 x 35 x 34 (22,000 compartments) and took 
less than 3 hours in the YeS-1045. 

6. Conclusion 


The paper describes a conservative method of unified com- 
putation of three-dimensional flow in the region of a complex 
shape. The accuracy suitable for practical purposes, the rela- 
tively modest computer run time, and the simple structure of the 
program complex render it a useful tool in the study of engineer- 
ing gas dynamic problems. 

The authors express deep gratitude to A. A. Samarskiy for 
unwavering attention to the work, as well as A. 0. Latsis for 
useful discussions of the elements of the program realization 
of the method. 
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